Meta-analysis of the Interoceptive Accuracy Scale

Study 1

Code
library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)

Data Preparation

Murphy (https://osf.io/3m5nh)

Code
df1a <- haven::read_sav("../data/Murphy2020/Study 1.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

df1b <- haven::read_sav("../data/Murphy2020/Study 6 IAS.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

Gaggero (2021) (https://osf.io/5x9sg)

Code
load("../data/Gaggero2021/data.rda")
df2 <- data |> 
  mutate(Gender = as.character(gender)) |> 
  rename(
    Age=age, Heart = `IAS 1`, Hungry = `IAS 2`, Breathing = `IAS 3`, Thirsty = `IAS 4`,
    Urinate = `IAS 5`, Defecate = `IAS 6`, Taste = `IAS 7`, Vomit = `IAS 8`, Sneeze = `IAS 9`,
    Cough = `IAS 10`, Temp = `IAS 11`, Sex_arousal = `IAS 12`, Wind = `IAS 13`, Burp = `IAS 14`,
    Muscles = `IAS 15`, Bruise = `IAS 16`, Pain = `IAS 17`, Blood_Sugar = `IAS 18`,
    Affective_touch = `IAS 19`, Tickle = `IAS 20`, Itch = `IAS 21`)
rm(data)

Campos (2022) - Study 1 (https://osf.io/j6ef3)

Code
df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )

Todd (2022) - https://osf.io/ms354/

Code
df4 <- haven::read_sav("../data/Todd2022/CompleteData_new.sav") |>
  rename_with(\(x) str_remove(x, "IAS_"), .cols = starts_with("IAS_")) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )

Arslanova (2022) - https://osf.io/mp3cy/

Code
df5 <- readxl::read_excel("../data/Arslanova2022/Murphy_Iacc_scale.xlsx")  |> 
  filter(str_detect(Question.Key, "IAC_")) |> 
  filter(str_detect(Question.Key, "-quantised")) |> 
  pivot_wider(names_from = Question.Key, values_from = Response) |>
  rename_with(\(x) str_remove(x, "IAC_"), .cols = starts_with("IAC_")) |> 
  rename_with(\(x) str_remove(x, "-quantised"), .cols = everything()) |> 
  rename(Sex_arousal = SexuallyAroused, Itch = Itchy, Sex_arousal = SexuallyAroused,
         Temp = HotCold, Tickle = Ticklish, Breathing= BreatheFast, 
         Affective_touch= PleasantAffectionate, Blood_Sugar= BloodSugar, 
         Muscles=TiredMuscles, Heart= FastHeart, Taste=Tastes) |> 
  select(-starts_with("Participant")) |> 
  mutate_all(as.numeric)

Brand (2022) - https://osf.io/xwz6g/

Code
df6 <- haven::read_sav("../data/Brand2022/LatentVariableApproach.sav") |> 
  select(-SERIAL) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

Brand (2023) - https://osf.io/3f2h6/

Code
df7a <- haven::read_sav("../data/Brand2023/Data_Giessen.sav") |> 
  select(-COUNTRY_OTHER) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(SEX == 1, "Male", ifelse(SEX == 2, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

df7b <- haven::read_sav("../data/Brand2023/Data_Mainz.sav") |> 
  select(-COUNTRY_OTHER, -EDUCATION_OTHER) |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )
  
  
df7c <- haven::read_sav("../data/Brand2023/Data_PotVie.sav") |> 
  select(-ID) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IA02_01, Hungry = IA02_02, Breathing = IA02_03, Thirsty = IA02_04,
    Urinate = IA02_05, Defecate = IA02_06, Taste = IA02_07, Vomit = IA02_08, Sneeze = IA02_09,
    Cough = IA02_10, Temp = IA02_11, Sex_arousal = IA02_12, Wind = IA02_13, Burp = IA02_14,
    Muscles = IA02_15, Bruise = IA02_16, Pain = IA02_17, Blood_Sugar = IA02_18,
    Affective_touch = IA02_19, Tickle = IA02_20, Itch = IA02_21
  )

Lin (2023) - https://osf.io/3eztd - Note: No tickle because same Chinese character

Code
df8a <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |>
  select(-sex) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |>
  rename(
    Age = age,
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH
  ) 


df8b <- haven::read_sav("../data/Lin2023/Study 2.sav") |>
  mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |>
  rename(
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH)

VonMohr (2023) - Study 3 (https://osf.io/7p9u5/)

Code
df9 <- read.csv("../data/VonMohr2023/DataSet_study3.csv")
df9 <- df9[complete.cases(select(df9, starts_with("IAS_"))), ]
df9 <- df9 |>
  mutate(Gender = as.character(ifelse(GENDER == "Man", "Male", ifelse(GENDER == "Woman", "Female", "Other")))) |>
  rename(
    Age=AGE, Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

Makowski

  • Note: First sample has some missing items
  • Note: Analog scales
Code
df10a <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
  rename(
    Gender = Sex, Heart = Item_IAS_Interoception_1, Hungry = Item_IAS_Interoception_2,
    Breathing = Item_IAS_Interoception_3, Thirsty = Item_IAS_Interoception_4,
    Temp = Item_IAS_Interoception_5, Sex_arousal = Item_IAS_Interoception_6,
    Urinate = Item_IAS_Elimination_1, Defecate = Item_IAS_Elimination_2,
    Vomit = Item_IAS_Elimination_3, Wind = Item_IAS_Expulsion_1,
    Burp = Item_IAS_Expulsion_2, Sneeze = Item_IAS_Expulsion_3,
    Muscles = Item_IAS_Nociception_1, Bruise = Item_IAS_Nociception_2,
    Pain = Item_IAS_Nociception_3, Affective_touch = Item_IAS_Skin_1,
    Tickle = Item_IAS_Skin_2, Itch = Item_IAS_Skin_3
  ) |>
  filter(!is.na(Urinate))

df10b <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
  rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

data <- list(df1a=df1a, df1b=df1b, df2=df2, df3=df3, df4=df4, df5=df5, df6=df6, df7a=df7a, df7b=df7b, df7c=df7c, df8a=df8a, df8b=df8b, df9=df9, df10a=df10a, df10b=df10b)
save(data, file = "../data/data.RData")

Participants

  • Sample 1a: Data from Murphy’s (2020) study 1, downloaded from OSF, included 451 participants (Mean age = 25.8, SD = 8.4, range: [18, 69]; Gender: 69.4% women, 29.5% men, 1.11% non-binary).
  • Sample 1b: Data from Murphy’s (2020) study 6, downloaded from OSF, included 375 participants (Mean age = 35.3, SD = 16.9, range: [18, 91]; Gender: 70.1% women, 28.5% men, 1.33% non-binary).
  • TODO: ADD INFO

Total N = 31317.

Descriptive

Distribution

The items with the differing distribution pattern (with a lower mode around 2/5) are “Affective Touch”, “Blood Sugar” and “Bruise”.

Code
vars <- c(
  "Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
  "Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch"
)

dens1a <- estimate_density(select(df1a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1a")
dens1b <- estimate_density(select(df1b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1b")
dens2 <- estimate_density(select(df2, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 6")
dens7a <- estimate_density(select(df7a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7a")
dens7b <- estimate_density(select(df7b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7b")
dens7c <- estimate_density(select(df7c, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7c")
dens8a <- estimate_density(select(df8a, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8a")
dens8b <- estimate_density(select(df8b, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8b")
dens9 <- estimate_density(select(df9, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 9")
dens10a <- estimate_density(select(df10a, all_of(setdiff(vars, c("Taste", "Cough", "Blood_Sugar")))), method = "kernSmooth") |>
  mutate(Sample = "Sample 10a",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10b <- estimate_density(select(df10b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 10b",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))


p1 <- rbind(dens1a, dens1b, dens2, dens3, dens4, dens5, dens6, dens7a, dens7b, dens7c, dens8a, dens8b, dens9, dens10a, dens10b) |>
  mutate(Sample = fct_relevel(Sample, "Sample 1a", "Sample 1b", "Sample 2", "Sample 3", "Sample 4", "Sample 5", "Sample 6", "Sample 7a", "Sample 7b", "Sample 7c", "Sample 8a", "Sample 8b", "Sample 9", "Sample 10a", "Sample 10b")) |>
  # mutate(Dashed = ifelse(Parameter %in% c("Affective_touch", "Blood_Sugar", "Bruise"), TRUE, FALSE)) |> 
  ggplot(aes(x = x, y = y, color = Parameter)) +
  geom_line() +
  scale_color_material_d() +
  scale_linetype_manual(values = c("solid", "dashed"), guide="none") +
  theme_minimal() +
  theme(axis.text.y = element_blank(),
        plot.title = element_text(face="bold")) +
  labs(title = "Distribution of responses for all items across various datasets", x = "Response", y = "Density", color = "Item") +
  facet_wrap(~Sample, scales = "free_y", nrow = 5)

ggsave("figures/Figure1.png", p1, width=10, height=10, dpi=200, bg="white")
p1

Code
data1a <- normalize(select(df1a, all_of(dens1a$Parameter)), range = c(1, 5))
data1b <- normalize(select(df1b, all_of(dens1b$Parameter)), range = c(1, 5))
data2 <- normalize(select(df2, all_of(dens2$Parameter)), range = c(1, 5))
data3 <- normalize(select(df3, all_of(dens3$Parameter)), range = c(1, 5))
data4 <- normalize(select(df4, all_of(dens4$Parameter)), range = c(1, 5))
data5 <- normalize(select(df5, all_of(dens5$Parameter)), range = c(1, 5))
data6 <- normalize(select(df6, all_of(dens6$Parameter)), range = c(1, 5))
data7a <- normalize(select(df7a, all_of(dens7a$Parameter)), range = c(1, 5))
data7b <- normalize(select(df7b, all_of(dens7b$Parameter)), range = c(1, 5))
data7c <- normalize(select(df7c, all_of(dens7c$Parameter)), range = c(1, 5))
data8a <- normalize(select(df8a, all_of(dens8a$Parameter)), range = c(1, 5))
data8b <- normalize(select(df8b, all_of(dens8b$Parameter)), range = c(1, 5))
data9 <- normalize(select(df9, all_of(dens9$Parameter)), range = c(1, 5))
data10a <- select(df10a, all_of(dens10a$Parameter))
data10b <- select(df10b, all_of(dens10b$Parameter))

data_all <- rbind(
  data1a, data1b, data2, data3, data4, data5, data6, data7a, data7b, data7c,
  mutate(data8a, Tickle = NA), mutate(data8b, Tickle = NA), data9,
  mutate(data10a, Taste = NA, Cough = NA, Blood_Sugar = NA), data10b
) 

Correlations

An overall positive intercorrelation pattern, with no clear structure emerging.

Code
make_correlation <- function(df) {
  correlation::correlation(df, redundant = TRUE) |>
    correlation::cor_sort() |>
    correlation::cor_lower() |>
    mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
    mutate(Parameter2 = fct_rev(Parameter2)) |>
    ggplot(aes(x = Parameter1, y = Parameter2)) +
    geom_tile(aes(fill = r), color = "white") +
    geom_text(aes(label = val), size = 3) +
    labs(title = "Correlation Matrix") +
    scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks = FALSE)) +
    theme_minimal() +
    theme(
      legend.title = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1)
    )
}

make_correlation(data_all)

Code
make_correlation(data1a)

Code
make_correlation(data1b)

Code
make_correlation(data2)

Code
make_correlation(data3)

Code
make_correlation(data4)

Code
make_correlation(data5)

Code
make_correlation(data6)

Code
make_correlation(data7a)

Code
make_correlation(data7b)

Code
make_correlation(data7c)

Code
make_correlation(data8a)

Code
make_correlation(data8b)

Code
make_correlation(data9)

Code
make_correlation(data10a)

Code
make_correlation(data10b)

Graph Analysis

Exploratory Graph Analysis (EGA) is a recently developed framework for psychometric assessment, that can be used to estimate the number of dimensions in questionnaire data using network estimation methods and community detection algorithms, and assess the stability of dimensions and items using bootstrapping. See https://r-ega.net/articles/ega.html for details.

Unique Variable Analysis (UVA)

Unique Variable Analysis (Christensen, Garrido, & Golino, 2023) uses the weighted topological overlap measure (Nowick et al., 2009) on an estimated network. Values greater than 0.25 are determined to have considerable local dependence (i.e., redundancy) that should be handled (variables with the highest maximum weighted topological overlap to all other variables (other than the one it is redundant with) should be removed).

Code
uva <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.368

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.293
 Urinate Defecate 0.270

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Sneeze     Cough 0.234
  Heart Breathing 0.225
 Hungry   Thirsty 0.218
Code
uva$keep_remove
$keep
[1] "Itch"

$remove
[1] "Tickle"
Code
uva <- EGAnet::UVA(data = data1a, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.286
   Wind   Burp 0.275

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.244
 Urinate Defecate 0.241
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data1b, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.350
 Sneeze  Cough 0.317
   Wind   Burp 0.309

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva$keep_remove
$keep
[1] "Cough" "Burp"  "Itch" 

$remove
[1] "Sneeze" "Wind"   "Tickle"
Code
uva <- EGAnet::UVA(data = data3, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.445
 Sneeze  Cough 0.318

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
   Wind   Burp 0.256

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i      node_j   wto
  Bruise Blood_Sugar 0.219
 Urinate    Defecate 0.217
   Heart   Breathing 0.208
Code
uva$keep_remove
$keep
[1] "Sneeze" "Itch"  

$remove
[1] "Cough"  "Tickle"
Code
uva <- EGAnet::UVA(data = data4, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.415
 Tickle      Itch 0.351

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
  Hungry  Thirsty 0.293
 Urinate Defecate 0.282
    Wind     Burp 0.275
  Sneeze    Cough 0.251

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i node_j   wto
 Muscles   Pain 0.205
Code
uva$keep_remove
$keep
[1] "Heart" "Itch" 

$remove
[1] "Breathing" "Tickle"   
Code
uva <- EGAnet::UVA(data = data5, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.251

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
 Bruise   Itch 0.241
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.416
  Sneeze    Cough 0.332
    Wind     Burp 0.314
  Tickle     Itch 0.303

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i  node_j   wto
 Hungry Thirsty 0.274

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.237
Code
uva$keep_remove
$keep
[1] "Urinate" "Sneeze"  "Wind"    "Itch"   

$remove
[1] "Defecate" "Cough"    "Burp"     "Tickle"  
Code
uva <- EGAnet::UVA(data = data7a, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Sneeze  Cough 0.434
 Tickle   Itch 0.376

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.295
 Urinate Defecate 0.294

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva$keep_remove
$keep
[1] "Sneeze" "Tickle"

$remove
[1] "Cough" "Itch" 
Code
uva <- EGAnet::UVA(data = data7b, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.317
  Tickle     Itch 0.311

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.281
 Sneeze     Cough 0.261

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.237
Code
uva$keep_remove
$keep
[1] "Defecate" "Itch"    

$remove
[1] "Urinate" "Tickle" 
Code
uva <- EGAnet::UVA(data = data7c, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.457
 Sneeze  Cough 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.280
 Urinate Defecate 0.254

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Hungry   Thirsty 0.242
  Vomit    Sneeze 0.224
  Heart Breathing 0.218
Code
uva$keep_remove
$keep
[1] "Cough" "Itch" 

$remove
[1] "Sneeze" "Tickle"
Code
uva <- EGAnet::UVA(data = data8a, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.237
   Heart Breathing 0.218
  Hungry   Thirsty 0.213
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data8b, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.277
   Heart Breathing 0.267

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.206
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data9, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
  Tickle     Itch 0.379
    Wind     Burp 0.373
 Urinate Defecate 0.341

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.285
 Sneeze     Cough 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j  wto
 Hungry Thirsty 0.24
Code
uva$keep_remove
$keep
[1] "Defecate" "Wind"     "Itch"    

$remove
[1] "Urinate" "Burp"    "Tickle" 
Code
uva <- EGAnet::UVA(data = data10a, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.297

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.209
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data10b, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva$keep_remove
NULL

Networks

Code
make_egas <- function(data) {
  egas <- list()
  for (model in c("glasso")) {  # , "TMFG"
    for (algo in c("walktrap", "louvain")) {
      for (type in c("ega", "ega.fit")) { # "" # "hierega"
        if (type == "ega.fit" & algo == "louvain") next # Too slow
        egas[[paste0(model, "_", algo, "_", type)]] <- EGAnet::bootEGA(
          data = data,
          seed = 123,
          model = model,
          algorithm = algo,
          EGA.type = type,
          type = "resampling",
          plot.itemStability = FALSE,
          verbose = FALSE
        )
      }
    }
  }

  p <- EGAnet::compare.EGA.plots(
    egas$glasso_walktrap_ega, 
    # egas$glasso_walktrap_ega.fit,
    egas$glasso_louvain_ega, 
    # egas$TMFG_louvain_ega,
    # egas$TMFG_walktrap_ega, 
    # egas$TMFG_walktrap_ega.fit,
    labels = c(
      "glasso_walktrap_ega", 
      # "glasso_walktrap_ega.fit",
      "glasso_louvain_ega" 
      # "TMFG_louvain_ega",
      # "TMFG_walktrap_ega", 
      # "TMFG_walktrap_ega.fit"
    ),
    rows = 3,
    plot.all = FALSE
  )$all
  list(egas = egas, p = p)
}


egas_all <- make_egas(data_all)
egas_all$p

Code
egas1a <- make_egas(data1a)
egas1a$p

Code
egas1b <- make_egas(data1b)
egas1b$p

Code
egas2 <- make_egas(data2)
egas2$p

Code
egas3 <- make_egas(data3)
egas3$p

Code
egas4 <- make_egas(data4)
egas4$p

Code
egas5 <- make_egas(data5)
egas5$p

Code
egas6 <- make_egas(data6)
egas6$p

Code
egas7a <- make_egas(data7a)
egas7a$p

Code
egas7b <- make_egas(data7b)
egas7b$p

Code
egas7c <- make_egas(data7c)
egas7c$p

Code
egas8a <- make_egas(data8a)
egas8a$p

Code
egas8b <- make_egas(data8b)
egas8b$p

Code
egas9 <- make_egas(data9)
egas9$p

Code
egas10a <- make_egas(data10a)
egas10a$p

Code
egas10b <- make_egas(data10b)
egas10b$p

Structure Stability

Code
patchwork::wrap_plots(lapply(egas_all$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas1a$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas1b$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas2$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas3$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas4$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas5$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas6$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas7a$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas7b$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas7c$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas8a$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas8b$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas9$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas10a$egas, plot), nrow = 3)

Code
patchwork::wrap_plots(lapply(egas10b$egas, plot), nrow = 3)

Scores

Code
ega <- egas_all$egas$glasso_louvain_ega$EGA
plot(ega)

Code
get_scores <- function(data, ega) {
  EGAnet::net.scores(data, ega)$scores$std.scores |>
    as.data.frame() |>
    setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))
}

ega_scores_1a <- get_scores(data1a, ega)
ega_scores_1b <- get_scores(data1b, ega)
ega_scores_2 <- get_scores(data2, ega)
ega_scores_3 <- get_scores(data3, ega)
ega_scores_4 <- get_scores(data4, ega)
ega_scores_5 <- get_scores(data5, ega)
ega_scores_6 <- get_scores(data6, ega)
ega_scores_7a <- get_scores(data7a, ega)
ega_scores_7b <- get_scores(data7b, ega)
ega_scores_7c <- get_scores(data7c, ega)
ega_scores_8a <- get_scores(data8a, ega)
ega_scores_8b <- get_scores(data8b, ega)
ega_scores_9 <- get_scores(data9, ega)
ega_scores_10a <- get_scores(data10a, ega)
ega_scores_10b <- get_scores(data10b, ega)

Factor Analysis

How Many Factors

Code
n <- parameters::n_factors(data_all, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
3 Scree (SE) Scree_SE
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data1a, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
3 Scree (SE) Scree_SE
3 BIC BIC
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data1b, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
3 CNG CNG
3 Scree (SE) Scree_SE
3 BIC BIC
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
5 BIC (adjusted) BIC
Code
n <- parameters::n_factors(data2, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
4 beta Multiple_regression
4 Scree (SE) Scree_SE
4 BIC BIC
5 Optimal coordinates Scree
5 Parallel analysis Scree
5 Kaiser criterion Scree
Code
n <- parameters::n_factors(data3, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
3 CNG CNG
3 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 Scree (SE) Scree_SE
4 BIC BIC
Code
n <- parameters::n_factors(data4, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 Bentler Bentler
3 CNG CNG
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 Scree (SE) Scree_SE
5 BIC BIC
Code
n <- parameters::n_factors(data5, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Optimal coordinates Scree
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
1 BIC BIC
3 CNG CNG
3 BIC (adjusted) BIC
4 Bartlett Barlett
4 beta Multiple_regression
4 Scree (SE) Scree_SE
5 Anderson Barlett
5 Lawley Barlett
Code
n <- parameters::n_factors(data6, n_max = 12)
Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully
Code
plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 Bentler Bentler
3 CNG CNG
3 Scree (SE) Scree_SE
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data7a, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
3 Scree (SE) Scree_SE
3 BIC BIC
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data7b, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
4 beta Multiple_regression
4 Scree (SE) Scree_SE
5 Bentler Bentler
5 Optimal coordinates Scree
5 Parallel analysis Scree
5 Kaiser criterion Scree
Code
n <- parameters::n_factors(data7c, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 Bentler Bentler
3 CNG CNG
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Scree (SE) Scree_SE
Code
n <- parameters::n_factors(data8a, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 Bentler Bentler
3 CNG CNG
4 beta Multiple_regression
4 Scree (SE) Scree_SE
5 Optimal coordinates Scree
5 Parallel analysis Scree
5 Kaiser criterion Scree
5 BIC BIC
Code
n <- parameters::n_factors(data8b, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
4 beta Multiple_regression
4 BIC BIC
5 Optimal coordinates Scree
5 Parallel analysis Scree
5 Kaiser criterion Scree
5 Scree (SE) Scree_SE
5 BIC (adjusted) BIC
Code
n <- parameters::n_factors(data9, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
3 Scree (SE) Scree_SE
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data10a, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
3 Scree (SE) Scree_SE
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 BIC BIC
4 BIC (adjusted) BIC
Code
n <- parameters::n_factors(data10b, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
3 CNG CNG
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 Scree (SE) Scree_SE
5 BIC (adjusted) BIC

EFA

Code
efa5 <- parameters::factor_analysis(data_all, n = 4, rotation = "oblimin", sort = TRUE)
Loading required namespace: GPArotation
Code
plot(efa5)

Code
display(efa5)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR2 MR4 Complexity Uniqueness
Burp 0.74 0.03 -0.02 -0.07 1.02 0.48
Cough 0.70 -0.02 0.02 0.06 1.02 0.47
Wind 0.67 0.02 -0.02 -2.95e-03 1.00 0.55
Sneeze 0.63 -0.06 0.04 0.14 1.12 0.52
Vomit 0.44 0.06 0.04 0.18 1.39 0.63
Temp 0.27 0.24 0.05 0.22 2.98 0.61
Sex_arousal 0.26 0.22 0.03 0.20 2.90 0.67
Taste 0.24 0.20 0.10 0.15 2.99 0.70
Breathing 0.05 0.61 -0.04 0.04 1.03 0.59
Hungry -0.11 0.54 -5.62e-03 0.19 1.34 0.64
Heart 0.08 0.49 0.03 -0.07 1.10 0.73
Thirsty -0.07 0.47 -7.04e-03 0.25 1.60 0.64
Pain 0.16 0.40 0.10 0.10 1.60 0.62
Muscles 0.30 0.36 0.07 0.03 2.03 0.59
Blood_Sugar 0.11 0.35 0.22 -0.18 2.47 0.76
Bruise 0.25 0.34 0.22 -0.19 3.30 0.65
Tickle -0.04 -0.04 0.82 0.04 1.02 0.37
Itch 0.02 6.24e-03 0.75 -0.04 1.01 0.43
Affective_touch 0.08 0.19 0.33 0.11 2.00 0.70
Urinate 0.04 0.10 0.04 0.65 1.06 0.45
Defecate 0.15 0.04 0.05 0.61 1.14 0.47

The 4 latent factors (oblimin rotation) accounted for 41.65% of the total variance of the original data (MR1 = 14.63%, MR3 = 11.59%, MR2 = 8.13%, MR4 = 7.30%).

Code
efa4 <- parameters::factor_analysis(data1a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR2 MR3 MR1 MR4 Complexity Uniqueness
Tickle 0.72 0.03 -0.10 0.11 1.09 0.47
Itch 0.67 0.14 -9.44e-03 -0.10 1.13 0.46
Affective_touch 0.49 -0.13 0.29 0.13 1.97 0.59
Blood_Sugar 0.29 0.21 0.18 -0.23 3.51 0.75
Burp -0.03 0.73 -0.03 -1.21e-03 1.00 0.50
Cough 7.56e-03 0.69 7.00e-03 -0.01 1.00 0.52
Sneeze 0.07 0.54 7.09e-03 0.20 1.30 0.55
Wind 0.08 0.46 0.07 0.16 1.37 0.62
Hungry -0.16 0.04 0.65 0.03 1.13 0.60
Thirsty -4.84e-03 0.01 0.49 0.11 1.11 0.70
Breathing 0.14 0.05 0.44 0.03 1.24 0.70
Heart 0.11 -0.07 0.40 0.08 1.32 0.78
Bruise 0.31 0.25 0.34 -0.25 3.74 0.59
Pain 0.29 0.01 0.32 0.20 2.67 0.61
Temp 0.19 0.13 0.29 0.16 2.91 0.67
Muscles 0.20 0.13 0.26 0.16 3.24 0.69
Sex_arousal 0.18 0.02 0.21 0.15 2.81 0.82
Defecate 0.10 0.11 2.77e-03 0.66 1.10 0.46
Urinate -0.11 0.13 0.24 0.52 1.65 0.54
Taste 0.23 0.06 0.17 0.30 2.60 0.69
Vomit 0.21 0.18 0.11 0.25 3.26 0.71

The 4 latent factors (oblimin rotation) accounted for 37.94% of the total variance of the original data (MR2 = 10.46%, MR3 = 10.17%, MR1 = 10.08%, MR4 = 7.23%).

Code
efa4 <- parameters::factor_analysis(data1b, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR2 MR3 MR4 Complexity Uniqueness
Urinate 0.69 -0.18 0.13 -0.04 1.22 0.54
Defecate 0.68 -0.13 0.04 0.06 1.10 0.55
Breathing 0.56 0.03 0.06 7.96e-03 1.03 0.63
Hungry 0.55 0.15 -7.15e-03 0.02 1.16 0.60
Thirsty 0.52 0.05 0.14 -0.01 1.16 0.62
Sex_arousal 0.50 0.14 -0.04 0.20 1.49 0.56
Temp 0.46 0.23 0.16 -0.01 1.77 0.52
Pain 0.46 0.43 -0.05 -0.03 2.03 0.49
Taste 0.40 0.21 -0.07 0.12 1.77 0.69
Heart 0.40 0.06 -0.04 0.15 1.36 0.76
Muscles 0.35 0.24 0.07 0.09 2.07 0.65
Vomit 0.32 0.06 0.23 0.11 2.20 0.68
Tickle 0.02 0.69 0.01 0.04 1.01 0.47
Itch -0.02 0.67 0.18 0.03 1.16 0.41
Bruise 0.07 0.44 0.15 0.06 1.32 0.65
Blood_Sugar -0.08 0.38 0.16 0.09 1.61 0.77
Affective_touch 0.38 0.38 -0.13 0.05 2.26 0.63
Sneeze 0.05 0.05 0.71 0.01 1.02 0.41
Cough 0.04 0.02 0.71 0.04 1.01 0.43
Wind 7.52e-03 -0.04 -5.18e-03 0.94 1.00 0.15
Burp -0.01 0.27 0.17 0.47 1.91 0.48

The 4 latent factors (oblimin rotation) accounted for 44.25% of the total variance of the original data (MR1 = 17.71%, MR2 = 11.72%, MR3 = 7.65%, MR4 = 7.17%).

CFA

We discarded the following items: - Tickle: not present in the same dataset and consistently flagged as redundant.

Ambiguous: - Temp - Vomit - Affective_touch - Sex_arousal - Taste

Lower-Level Factors

Code
m1 <- "
Interoception =~ Hungry + Thirsty + Urinate + Defecate + Itch + Bruise + Muscles + Pain + Breathing + Heart + Cough + Sneeze + Wind + Burp
"

m4 <- "
Sustenance =~ Hungry + Thirsty + Urinate + Defecate + Muscles + Pain
Skin =~ Itch + Bruise
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"


m5 <- "
Sustenance =~ Hungry + Thirsty + Pain
UrinateDefecate =~ Urinate + Defecate
Skin =~ Itch + Bruise + Muscles
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"


m6 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"

m7 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
CoughSneeze =~ Cough + Sneeze
WindBurp =~ Wind + Burp
"

is_converged <- function(m) {
  ev <- eigen(lavaan::lavTech(m, "cov.lv")[[1]], symmetric=TRUE, only.values=TRUE)$values
  if (any(ev < -1 * .Machine$double.eps^(3/4))) {
    "Not converged"
  } else {
    ""
  }
}

compare_mods <- function(data) {
  m1 <- lavaan::cfa(m1, data = data)
  m4 <- lavaan::cfa(m4, data = data)
  m5 <- lavaan::cfa(m5, data = data)
  m6 <- lavaan::cfa(m6, data = data)
  m7 <- lavaan::cfa(m7, data = data)
  
  anova(m1, m4, m5, m6, m7) |>
    parameters::parameters() |>
    mutate(AIC = format(round(AIC, 0), nsmall = 0),
           BIC = format(round(BIC, 0), nsmall = 0),
           Chi2 = format(round(Chi2, 2), nsmall = 2),
           Converged = sapply(list(m1, m4, m5, m6, m7), is_converged))
}

display(compare_mods(data_all))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -168983 -168574 2080.93 56
m6 -164335 -163977 6740.02 62 < .001
m5 -160047 -159731 11038.16 67 < .001
m4 -154973 -154690 16120.58 71 < .001
m1 -141854 -141621 29251.50 77 < .001
Code
m7_all <- lavaan::cfa(m7, data = data_all)
parameters::parameters(m7_all, standardize=TRUE) |>
  display()
# Loading
Link Coefficient SE 95% CI z p
HungryThirsty =~ Hungry 0.68 4.79e-03 (0.67, 0.69) 142.07 < .001
HungryThirsty =~ Thirsty 0.71 4.77e-03 (0.70, 0.72) 148.26 < .001
UrinateDefecate =~ Urinate 0.77 3.94e-03 (0.76, 0.78) 195.43 < .001
UrinateDefecate =~ Defecate 0.77 3.93e-03 (0.77, 0.78) 196.53 < .001
ItchBruise =~ Itch 0.50 5.73e-03 (0.49, 0.51) 87.60 < .001
ItchBruise =~ Bruise 0.67 5.99e-03 (0.66, 0.68) 112.01 < .001
MusclesPain =~ Muscles 0.71 4.15e-03 (0.70, 0.72) 171.30 < .001
MusclesPain =~ Pain 0.68 4.22e-03 (0.68, 0.69) 161.84 < .001
HeartBreathing =~ Breathing 0.79 5.19e-03 (0.78, 0.80) 151.51 < .001
HeartBreathing =~ Heart 0.59 5.13e-03 (0.58, 0.60) 114.83 < .001
CoughSneeze =~ Cough 0.81 3.49e-03 (0.80, 0.81) 230.73 < .001
CoughSneeze =~ Sneeze 0.75 3.65e-03 (0.74, 0.76) 205.71 < .001
WindBurp =~ Wind 0.77 3.58e-03 (0.76, 0.78) 215.22 < .001
WindBurp =~ Burp 0.82 3.44e-03 (0.81, 0.83) 238.24 < .001
# Correlation
Link Coefficient SE 95% CI z p
HungryThirsty ~~ UrinateDefecate 0.65 6.26e-03 (0.63, 0.66) 103.09 < .001
HungryThirsty ~~ ItchBruise 0.45 9.03e-03 (0.43, 0.46) 49.41 < .001
HungryThirsty ~~ MusclesPain 0.63 7.04e-03 (0.61, 0.64) 89.04 < .001
HungryThirsty ~~ HeartBreathing 0.64 7.01e-03 (0.63, 0.66) 91.55 < .001
HungryThirsty ~~ CoughSneeze 0.45 7.14e-03 (0.44, 0.47) 63.70 < .001
HungryThirsty ~~ WindBurp 0.43 7.12e-03 (0.42, 0.44) 60.47 < .001
UrinateDefecate ~~ ItchBruise 0.43 8.38e-03 (0.42, 0.45) 51.81 < .001
UrinateDefecate ~~ MusclesPain 0.64 6.23e-03 (0.63, 0.66) 103.29 < .001
UrinateDefecate ~~ HeartBreathing 0.53 6.81e-03 (0.52, 0.55) 78.24 < .001
UrinateDefecate ~~ CoughSneeze 0.56 6.03e-03 (0.54, 0.57) 92.16 < .001
UrinateDefecate ~~ WindBurp 0.50 6.19e-03 (0.49, 0.52) 81.27 < .001
ItchBruise ~~ MusclesPain 0.80 7.96e-03 (0.79, 0.82) 100.88 < .001
ItchBruise ~~ HeartBreathing 0.51 8.74e-03 (0.49, 0.53) 58.21 < .001
ItchBruise ~~ CoughSneeze 0.64 7.69e-03 (0.62, 0.65) 82.74 < .001
ItchBruise ~~ WindBurp 0.64 7.55e-03 (0.62, 0.65) 84.63 < .001
MusclesPain ~~ HeartBreathing 0.64 6.99e-03 (0.62, 0.65) 91.42 < .001
MusclesPain ~~ CoughSneeze 0.66 6.06e-03 (0.65, 0.68) 109.42 < .001
MusclesPain ~~ WindBurp 0.63 6.10e-03 (0.62, 0.65) 103.72 < .001
HeartBreathing ~~ CoughSneeze 0.52 6.81e-03 (0.50, 0.53) 75.70 < .001
HeartBreathing ~~ WindBurp 0.45 6.95e-03 (0.44, 0.47) 65.21 < .001
CoughSneeze ~~ WindBurp 0.72 4.84e-03 (0.71, 0.73) 149.24 < .001
Code
display(compare_mods(data1a))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -1951 -1749 88.22 56
m6 -1925 -1748 126.29 62 < .001
m5 -1860 -1703 201.38 67 < .001 Not converged
m4 -1834 -1694 234.82 71 < .001
m1 -1661 -1546 420.24 77 < .001

Conclusion: The 7th factor solution was preferred in most datasets, including by indices that penalize increased number of parameters (such as BIC).

Higher Order Factors

Code
m7h1 <- paste0(m7, "
              Interoception =~ HungryThirsty + UrinateDefecate + ItchBruise + MusclesPain + HeartBreathing + CoughSneeze + WindBurp
              ")

m7h2 <- paste0(m7, "
              Expulsion =~ CoughSneeze + WindBurp
              Interoception =~ HungryThirsty + UrinateDefecate + HeartBreathing
              ")

m7h3 <- paste0(m7, "
              Valenced =~ MusclesPain + HeartBreathing + ItchBruise
              Expulsion =~ CoughSneeze + WindBurp
              Homeostasis =~ HungryThirsty + UrinateDefecate
              ")

compare_mods_h <- function(data) {
  m7 <- lavaan::cfa(m7, data = data)
  m7h2 <- lavaan::cfa(m7h2, data = data)
  m7h1 <- lavaan::cfa(m7h1, data = data)
  m7h3 <- lavaan::cfa(m7h3, data = data)
  
  anova(m7, m7h1, m7h2, m7h3) |>
    parameters::parameters() |>
    mutate(AIC = format(round(AIC, 0), nsmall = 0),
           BIC = format(round(BIC, 0), nsmall = 0),
           Chi2 = format(round(Chi2, 2), nsmall = 2),
           Converged = sapply(list(m7, m7h1, m7h2, m7h3), is_converged))
}

display(compare_mods_h(data_all))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -168983 -168574 2080.93 56
m7h2 -168224 -167899 2859.83 66 < .001
m7h3 -166932 -166615 4153.51 67 < .001
m7h1 -164501 -164210 6590.01 70 < .001

No evidence for higher order factors.

Code
display(compare_mods_h(data1a))
Fixed Effects
Parameter AIC BIC Chi2 df p Converged
m7 -1951 -1749 88.22 56
m7h3 -1904 -1747 157.51 67 < .001
m7h1 -1873 -1729 193.71 70 < .001

Model Performance

Code
rbind(
  performance::performance(m7_all, metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "All"),
  performance::performance(update(m7_all, data = data1a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 1a"),
  performance::performance(update(m7_all, data = data1b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 1b"),
  performance::performance(update(m7_all, data = data2), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 2"),
  performance::performance(update(m7_all, data = data3), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 3"),
  performance::performance(update(m7_all, data = data4), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 4"),
  performance::performance(update(m7_all, data = data5), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 5"),
  performance::performance(update(m7_all, data = data6), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 6"),
  performance::performance(update(m7_all, data = data7a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7a"),
  performance::performance(update(m7_all, data = data7b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7b"),
  performance::performance(update(m7_all, data = data7c), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7c"),
  performance::performance(update(m7_all, data = data8a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 8a"),
  performance::performance(update(m7_all, data = data8b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 8b"),
  performance::performance(update(m7_all, data = data9), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 9"),
  performance::performance(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data10a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 10a"),
  performance::performance(update(m7_all, data = data10b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 10b")
) |>
  display()
Chi2 RMSEA CFI SRMR Sample
2080.93 0.03 0.98 0.02 All
88.22 0.04 0.98 0.03 Sample 1a
66.50 0.02 0.99 0.03 Sample 1b
123.77 0.04 0.97 0.03 Sample 2
114.03 0.04 0.98 0.03 Sample 3
81.35 0.02 0.99 0.02 Sample 4
83.52 0.05 0.92 0.06 Sample 5
152.88 0.05 0.97 0.03 Sample 6
106.90 0.04 0.97 0.04 Sample 7a
188.21 0.03 0.98 0.02 Sample 7b
122.56 0.04 0.97 0.03 Sample 7c
115.89 0.03 0.98 0.02 Sample 8a
83.61 0.03 0.98 0.03 Sample 8b
2003.31 0.04 0.98 0.02 Sample 9
82.59 0.05 0.97 0.03 Sample 10a
91.47 0.05 0.97 0.04 Sample 10b

Scores

Note the usage of descriptive factor names relating directly to the items to avoid abstraction.

Code
# Add empirical variables
add_empirical <- function(data, df) {
  x <- data |>
    mutate(
      HungryThirsty = (df$Hungry + df$Thirsty) / 2,
      UrinateDefecate = (df$Urinate + df$Defecate) / 2,
      MusclesPain = (df$Muscles + df$Pain) / 2,
      HeartBreathing = (df$Heart + df$Breathing) / 2,
      ItchBruise = (df$Itch + df$Bruise) / 2,
      WindBurp = (df$Wind + df$Burp) / 2)
  if("Cough" %in% names(df)) {
    mutate(x, CoughSneeze = (df$Cough + df$Sneeze) / 2)
  } else {
    mutate(x, CoughSneeze = df$Sneeze)
  }
}

# Deal with missing data
temp <- predict(m7_all, newdata=data7c)
scores7c <- data.frame(matrix(ncol = 7, nrow = nrow(ega_scores_7c))) 
scores7c[!is.na(ega_scores_7c$EGA1), ] <- temp
names(scores7c) <- names(as.data.frame(temp))

scores <- list(
  sample1a = cbind(ega_scores_1a, data_addprefix(as.data.frame(predict(m7_all, newdata=data1a)), "CFA_")) |>
    add_empirical(data1a) |> 
    mutate(Sample = "Sample 1a"),
  sample1b = cbind(ega_scores_1b, data_addprefix(as.data.frame(predict(m7_all, newdata=data1b)), "CFA_")) |>
    add_empirical(data1b) |> 
    mutate(Sample = "Sample 1b"),
  sample2 = cbind(ega_scores_2, data_addprefix(as.data.frame(predict(m7_all, newdata=data2)), "CFA_")) |>
    add_empirical(data2) |> 
    mutate(Sample = "Sample 2"),
  sample3 = cbind(ega_scores_3, data_addprefix(as.data.frame(predict(m7_all, newdata=data3)), "CFA_")) |>
    add_empirical(data3) |> 
    mutate(Sample = "Sample 3"),
  sample4 = cbind(ega_scores_4, data_addprefix(as.data.frame(predict(m7_all, newdata=data4)), "CFA_")) |>
    add_empirical(data4) |> 
    mutate(Sample = "Sample 4"),
  sample5 = cbind(ega_scores_5, data_addprefix(as.data.frame(predict(m7_all, newdata=data5)), "CFA_")) |>
    add_empirical(data5) |> 
    mutate(Sample = "Sample 5"),
  sample6 = cbind(ega_scores_6, data_addprefix(as.data.frame(predict(m7_all, newdata=data6)), "CFA_")) |>
    add_empirical(data6) |> 
    mutate(Sample = "Sample 6"),
  sample7a = cbind(ega_scores_7a, data_addprefix(as.data.frame(predict(m7_all, newdata=data7a)), "CFA_")) |>
    add_empirical(data7a) |> 
    mutate(Sample = "Sample 7a"),
  sample7b = cbind(ega_scores_7b, data_addprefix(as.data.frame(predict(m7_all, newdata=data7b)), "CFA_")) |>
    add_empirical(data7b) |> 
    mutate(Sample = "Sample 7b"),
  sample7c = cbind(ega_scores_7c, data_addprefix(scores7c, "CFA_")) |> 
    add_empirical(data7c) |> 
    mutate(Sample = "Sample 7c"),
  sample8a = cbind(ega_scores_8a, data_addprefix(as.data.frame(predict(m7_all, newdata=data8a)), "CFA_")) |>
    add_empirical(data8a) |> 
    mutate(Sample = "Sample 8a"),
  sample8b = cbind(ega_scores_8b, data_addprefix(as.data.frame(predict(m7_all, newdata=data8b)), "CFA_")) |>
    add_empirical(data8b) |> 
    mutate(Sample = "Sample 8b"),
  sample9 = cbind(ega_scores_9, data_addprefix(as.data.frame(predict(m7_all, newdata=data9)), "CFA_")) |>
    add_empirical(data9) |> 
    mutate(Sample = "Sample 9"),
  sample10a = cbind(ega_scores_10a, data_addprefix(as.data.frame(predict(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data_all), newdata=data10a)), "CFA_")) |>
    add_empirical(data10a) |> 
    mutate(Sample = "Sample 10a", CFA_CoughSneeze = NA),
  sample10b = cbind(ega_scores_10b, data_addprefix(as.data.frame(predict(m7_all, newdata=data10b)), "CFA_")) |>
    add_empirical(data10b) |> 
    mutate(Sample = "Sample 10b")
)

save(scores, file = "../data/scores.RData")